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ABSTRACT 

It is shown that gravitational lensing by point masses arranged in an infinitely extended regular 
lattice can be studied analytically using the Weierstrass functions. In particular, we draw the 
critical curves and the caustic networks for the lenses arranged in regular-polygonal - square, 
equilateral triangle, regular hexagon - grids. From this, the mean number of positive parity 
images as a function of the average optical depth is derived and compared to the case of 
the infinitely extended field of randomly distributed lenses. We find that the high degree of 
the symmetry in the lattice arrangement leads to a significant bias towards canceling of the 
shear caused by the neighboring lenses on a given lens position and lensing behaviour that is 
qualitatively distinct from the random star field. We also discuss some possible connections 
to more realistic lensing scenarios. 



1 INTRODUCTION 

Beyond the Galactic and Local-Group (in particular, Magellanic 
Clouds and M31) microlensing experiments (Alcock et al. 2000; 
Afonso et al. 2003; Calchi Novati et al. 2005; de Jong et al. 2006), 
the main interest on the point-mass lenses lies in the effect of 
the individual stars in a lensing galaxy (Chang & Refsdal 1979, 
1984; Nityananda & Ostriker 1984; Irwin et al. 1989; Witt et al. 
1995; Schechter & Wambsganss 2002; Keeton et al. 2006). For 
these cases, each point-mass lens cannot be treated individually 
and their collective effects usually differ drastically from that of 
the simple linear superposition of them. The usual approach to 
this "high-optical-depth microlensing" problem is "inverse ray trac- 
ing" (Kayser et al. 1986; Schneider & Weiss 1987; Wambsganss et 
al. 1990, 1992), that is, examining statistical properties of lensing 
observables using Monte-Carlo realizations of the lensing system. 
However, as the number of the individual point-mass lenses to re- 
produce the realistic scenario well approaches the number of stars 
in the galaxy, this becomes very expensive rather quickly in terms 
of the required resources. Another complementary approach to the 
problem is from the random walk and the probability theory to- 
gether with the thermodynamic approximation (that is, effectively 
considering the limit at the infinite number of the stars), which has 
been quite successful in certain well-defined problems (Nityananda 
& Ostriker 1984; Schneider 1987; Schneider et al. 1992). 

In the following, we consider a different approach to the prob- 
lem: seeking complete analytic solutions. In many complex sys- 
tems, imposing high level of symmetry can reduce the problem to 
a simpler one, which may be otherwise too complicated to analyze. 
In a similar spirit, in this paper we consider gravitational lensing 
by equal point masses that are arranged in an infinitely extended 
regular lattice. After we derive the appropriate lens equation in 
Sect. 2, they are examined in Sect. 3 to find the corresponding crit- 
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ical curves and the caustic networks. In Sect. 4, we study the mean 
number of positive parity images, which is compared to the case 
of lensing by an infinite random star field. In Sect. 5, we consider 
the effect of the external potential by mean of adding external shear 
(and convergence), and discuss some possible connections to more 
realistic lensing scenarios. 

We argue that, while the system considered here may be some- 
what artificial in its construct and rather abstract in its nature, the 
study presented in this paper can lead us to some insight on mi- 
crolensing at high optical depth. 



2 LENS EQUATION 

In most astronomical situations, gravitational lensing is described 
by the lens equation y = x - Vi^, which is derived from the 
lowest-order nontrivial approximation to the path of light propa- 
gation (Schneider et al. 1992; Petters et al. 2001; Kochanek et al. 
2005). Here, y and x are the (2-dimensional) vectors representing 
the lines of sight toward, respectively, the angular positions of the 
source under the absence of lensing and the lensed image. The lens- 
ing potential if/ is the gravitational potential (equivalently, the grav- 
itational 'Shapiro' time delay) integrated along the line of sight up 
to a scale constant (Blandford & Narayan 1986). 

In particular, gravitational lensing by a point mass is described 
by the potential if/(x) = #| In \x - l\ if the point-mass lens is located 
along the line of sight indicated by the vector / (Paczyfiski 1986). 
Here, the scale constant 6 E is known as the 'Einstein (ring) radius' 
and determined by the mass of the lens (6\ oc M) and the distances 
among the source, the lens, and the observer. With this potential, 
we find the lens equation for point-mass lensing, 



x-l 

\x - l\ 2 



(1) 



where the unit of angular measurements has been chosen such that 
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6 E = 1. The equation is generalized to the case of multiple point- 
mass lenses with a smoothly-varying 'large scale' potential (Young 
1981) such that 



y = x ~Tj 



m ; x-li 



(2) 



Here, the unit of angular measurements is now given by 6% corre- 
sponding to a fiducial mass scale M and subsequently the 'external' 
lensing potential (fr ext should be rescaled appropriately. In addition, 
the masses of individual lenses m, enter into the equation as the 
ratio to the fiducial scale. 

Next, let us think of the case that the equal point-mass lenses 
are located at the lattice points on an infinitely extended square 
grid. In addition, we ignore the effect of the external potential for 
the moment. Utilizing the complex-number notation (Bourassa & 
Kantowski 1975; Subramanian et al. 1985; Witt 1990; An 2005), 
the lens equation (2) then can be written as 



n = z - f(z) ; 



i 



(3) 



2(m + nijw 



where rj, z, and A are the complexified variables corresponding to y, 
x, and /, respectively, and the over-bar notation is used to indicate 
complex conjugation. We further assume that the lens locations are 
given by A = 2{m + n\)u) where m and n are any pair of integers 
(both running from the negative infinity to the positive infinity) and 
oj is the half distance between adjacent grid points. Without any 
loss of generality, w can be restricted to be a positive real. Strictly 
speaking, the 'deflection function' f(z) given formally in terms of 
the infinite sum in equation (3) is not well-defined as it is not con- 
vergent. This is in fact a natural consequence of having an infinite 
total lensing mass, which is unphysical. However, what is impor- 
tant in our understanding of the lensing is not the absolute amount 
of the deflection but the difference of the deflections between neigh- 
bouring lines of sights. That is to say, if we add any constant to the 
deflection function in the lens equation, the resulting lens equation 
is reduced to the equivalent one without the constant by introduc- 
ing the "offset" to the source position to cancel the constant. Two 
equations then are observationally equivalent because there is no a 
priori information regarding the source position in the absence of 
the lens. Similarly, despite the seeming global failure of equation 
(3), we can still proceed by focusing on the local properties of the 
lensing described by equation (3). 

For this, we first write down the (formal) second-order com- 
plex derivative of the deflection function; 



m,n=—oo 

■ p'(z\a>,ico) 



[z - 2(m + n\)(jj] 3 
p'(z;g 2 ,0) 



(4) 



where 



g2 



-(?)' 



1/4 
167T 4 



(5) 



We find that this is related to the definition of the special function 
known as the Weierstrass elliptic function 1 (Abramowitz & Stegun 



Throughout this paper, the arguments of the Weierstrass funetions followed by a vertical bar - | - 
denote the half-periods whereas the elliptic invariants are indicated by those followed by a semicolon 
(;). For details, see the listed references or Appendix. Whenever it is appropriate, these arguments 
may be suppressed, provided that there is little danger of confusion. 



1972, or see also Appendix). Here, T x = T(x) is the gamma func- 
tion, p(z) is the Weierstrass elliptic function and p'(z) is its complex 
derivative. We also use an abbreviated notation for (the power to) 
the gamma function such that [r(x)]" = T" in order to avoid pro- 
liferation of parentheses and brackets. In addition, a>o ~ 1.854 is 
the real half-period of the lemniscatic case of p-function (see Ap- 
pendix) and K t = 7t/(2ai) 2 is the optical depth (Vietri & Ostriker 
1983; Paczynski 1986) of the lenses. If we consider the square 
cell given by the centres of each lensing lattice, it is straightfor- 
ward to see that there exists a lens corresponding to each cell with 
side length of 2a>. Hence, the mean surface density of point masses 
is in fact given by <x = (2o))' 2 , and therefore, the optical depth 
by k+ = ncr = n/(2a>) 2 . Moreover, the 90° -rotational symme- 
try of the system implies that the total shear due to all the other 
lenses on any given lens cancels out to be zero [mathematically 
this indicates the constant term in the Laurent-series expansion 
of f'(z) is null]. Consequently, we find that f'(z) = -p(z;g2,0) 
and f(z) = £(z;g2,0) + C. Here, ((z) is the Weierstrass zeta func- 
tion (Abramowitz & Stegun 1972) and C is an arbitrary complex 
constant. If we choose the 'centre' of the system to coincide with 
the lens at the coordinate origin, then we have C = 0. We note 
that, while the constant C cannot be independently determined, its 
choice is physically inconsequential. The situation is analogous to 
lensing by an infinitely extended mass screen. While a naive ex- 
pectation from the symmetry appears to indicate null deflections 
along every line of sight, detailed physical consideration leads to 
the uniform focusing with respect to an arbitrary choice of the cen- 
tre. In the current situation, the discrete translation symmetry im- 
plies the arbitrariness in the choice of the centre and consequently 
that of C, which results in the infinite sum in equation (3) not being 
well-determined. However, the constant-deflection term in the lens 
equation again only leads to a constant offset between the source 
and the image/lens plane, which has no observable consequence 
and therefore can be ignored. 

Some results are immediate from the resulting lens equation 



(6) 



(V) 



t) = z- f(z;g 2 ,0). 

For instance, the property of ((z) such that 

f [z + 2(m + m)w; g 2 , 0] = ( (z; g 2 , 0) + — (m - ni) 

2co 

for arbitrary integers m and n indicates that the vertices z, z + 2oj, 
z + 2(1 + i)aj, and z + 2iaj of a square in the lens plane map to the 
points on the source plane: i] = z- <T(z), r\ + N, r\ + (1 + i)N, and 
r\ + iN where K = 2ai - n/(2a>), and therefore the mean inverse 
magnification is given by [X/(2w)] 2 = [1 - n/(2a>) 2 ] 2 . The result 
confirms the expectation that the sufficiently large beam of light 
would see the system as if it were a uniform screen of mass with a 
convergence (or the optical depth) k = ncr = n/(2w) 2 . Along the 
real line, equation (6) maps each segment between two adjacent 
lenses onto the whole real line, and thus it is easily argued that 
there are infinite number of images. However, unless k+ = 1, all 
but a finite number of images have negative parity. 



2.1 lenses on triangular and hexagonal grid 

Analogous to the preceding case, we can also set up the lens 
equation for the equal point-mass lenses located at the lattice 
points on an infinitely extended equilateral-triangular grid. The 
60° -rotational symmetry of the system again implies that the to- 
tal shear on any given lens caused by the remaining set of lenses 
also cancels out. This fact, combined with the general definition 
of p-function indicates that the lens equation can again be written 
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Figure 1. Contour plots of the equipotential lines of the lensing potential for the lattice lenses. Top Left: square lattice. Top Right: triangular lattice. Bottom 
Left: hexagonal lattice with equation (14). Bottom Right: hexagonal lattice with a lens being at the center. 



down using the Weierstrass functions. If the lenses are placed over 
the grid given by A = 2ai{m + ne ni/3 ) where m and n are integers, 
the corresponding lens equation is given by 



r\ = z - f(z) ; f(z) = ({z\oj, e m/3 co) = ffc 0, g 3 ) 



where 



=m'= 



V5r 



1/3 



8tt 3 



(8) 



(9) 



and ll>2 ~ 1.530 is the real half-periods of the equiharmonic case of 
^-function (see Appendix) whereas K t = n/(2 V3oi 2 ) is the corre- 
sponding optical depth. Likewise, the quasiperiodicity of f (z; 0, gi), 



convergence is K t = n/(2 y/3w 2 ) (and the mean surface density <x 
of the point masses on the triangular grid is given by cr _1 =2 y/Sco 2 
where 2co is the side length of the unit triangular cell). 

We note that the the preceding two cases of the lens arrange- 
ment correspond to two of the three possible regular tessellations of 
the two-dimensional plane. It may be of some interest to consider 
the regular lens placement corresponding to the remaining regu- 
lar tessellation - the hexagonal or honeycomb tiling. It turns out 
that the regular hexagonal grid case is closely related to the equi- 
lateral triangular grid. 2 By removing every third of the lens that 
forms a (30° -rotated) larger triangular grid of side length of 2 V3w 
from the base triangular grid with side length of 2u>, the remaining 
lenses form a hexagonal grid (or the vertices of honeycomb cells). 



t[z + 2(m + ne 7Z " 3 )w;0,g 3 ] = <r(z;0,g 3 ) + (m + mT m/3 ) (10) 

where m and n are arbitrary integers, indicates that the vertex points 
z,z + 2oj, and z + 2e m/3 a> of the equilateral triangle in the lens plane 
map onto the source plane points: r/ = z- f (z), r/ + 2 and r\ + 3e m/3 
where 3 = 2a> - n/( V3w) so that the mean inverse magnification 
is again given by (1 - k*) 2 where the optical depth or the mean 



Note that they are in fact so-called dual tiling of each other. That is, if we consider a grid with 
vertex points at the centres of each triangular cell of the equilateral triangular grid, then the resulting 
lattice forms a regular-hexagonal (honeycomb) grid and also vice versa. 
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Figure 2. Critical curves (dotted lines) and caustic networks (solid lines) 
for lenses on a square lattice with a) = 1.4 [k* = n/(2co) 2 a 0.401]. 



Figure 3. Same as Fig. 2 except for a) = 1.2 [«•* = n/(2u>) ~ 0.545]. 



Consequently, the corresponding lens equation may be written as 



77 = z- 6feg 3 ); 

b(z;g 3 ) = <T(z|w,e m/3 w) - £(z| V3e m/ V V3iw) 
= f(«0,&)-«(r,0,-g) 



(11) 



where #3 is related to a> through equation (9) provided that 2u> is 
still the side length of the base triangular grid (and also that of the 
unit honeycomb cell). However, since every third of lens has been 
removed from the base triangular cell, the optical depth is reduced 
to k+ = 7t/(3 Via) 2 ), and thus 



gi 



3 3/2 r 



L/3 



2 4 7T 3 



(12) 



In addition, b'(z;g 3 ) = -b(z;g 3 ) (the primed symbol indicates the 
complex differentiation, i.e., the derivative with respect to z), The 
function h(z; g 3 ) is studied in Appendix A2. 1 . It is again straightfor- 
ward to establish that the mean inverse magnification is (1 - k*) 2 - 
Unlike the preceding two cases, the system described by lens equa- 
tion (11) does not have the lens at the centre, but the centre of 
the system corresponds to the location on the lens plane such that 
6'(0) = b(0) = 0. However, we note that this is deliberately chosen 
for mathematical simplicity, and does not imply any intrinsic differ- 
ence of the honeycomb grid from the square or triangular grid. In 
particular, we note the 120° -rotational symmetry of the system with 
respect to any lens location, which actually implies that any lens 
on the honeycomb grid experiences zero shear from the remaining 
point masses. In other words, for all three cases of the point masses 
on the regular grid, the total shear from all the adjacent lenses ex- 
actly cancels out at any lens position. 



2.2 lensing potential 

While most of the lensing properties of the lattice lens can be stud- 
ied using the lens equations, it is still useful to have an expression 



for the lensing potential for some purposes such as the time de- 
lay. Although the direct infinite sum of the individual logarithmic 
potential of the point-mass lens is divergent everywhere, one can 
get around this difficulty through the known antiderivative of the 
Weierstrass zeta function. First, we note the relation between the 
real potential and the complex deflection function (An 2005; An & 
Evans 2006); the complex lens equation is given by r\ = z - 2d,ift if 
the real potential is given by if/. Here, the operator 3= indicates the 
'Wirtinger derivative' (e.g., Schramm & Kayser 1995) with respect 
to z, that is, if if/ = if/(x, y) and z = x + yi, then 



dif/ dx dif/ dy 
dx dz dy dz 



1 dif/ i dif/ 
2~dx + 2 



Sv 



(13) 



because x = (z + z)/2 and y = i(z-z)/2. For the current scenario, we 
have /(z) = 29 ; (/r = 29 -1^ = 2d,if/ (note that if/ = if/ since if/ is real). 
Moreover, f(z) is a complex analytic function and so there exists a 
complex analytic function if/ c (z) such that f(z) = if/' c {z) = d-if/ c and 
d t ij/ c = 0. Then, d z (i(/ c +if/ c ) = d z if/ c + d z if/ c = d-if/ c + d=if/ c = /(z), and 
thus, 2if/ = if/ c + \j/ c , that is, if/(x, y) = %[if/ c {x + yi)] up to an additive 
constant (An 2005). Here, is the real part of a complex- valued 
function /. Consequently, we find that the real potential up to an 
additive constant for the regular lensing lattice is given by 



if/(x,y) 



In Hz; g 2 ,0)| 
o-fo 0,g 3 )\ 
cr(z; 0,g 3 ) 



o-(z;0,-g 3 /27) 



(14) 



where z = x + yi and cr(z;g2,g3) is the Weierstrass sigma function 
(Abramowitz & Stegun 1972), which can be defined as the anti-log- 
derivative of the Weierstrass zeta function, that is, (d/dz) In cr{z) = 
£(z). Here, we have also used the property of the complex logarithm 
function such that Dl[lng] = ln|g| where g is any complex- valued 
function. Note that an additional linear term may be required to 
be added in the expression for the potential if the centre of the sys- 
tem is chosen differently from those of the lens equations discussed 
earlier. Fig. 1 shows contour plots for the equipotential lines for the 
potential given in equations (14). 
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Figure 4. Caustic networks for a square lattice lens with ai - 
a) = 0.73 («* » 1.474; bottom right). 



1 (Ki, a 0.785; top left), at = 0.8 («* a 1.227; top right), ai = 0.77 (k* ss 1.325; bottom left), and 



3 CRITICAL CURVES AND CAUSTICS 

As usual, the Jacobian determinant of equation (6), (8) or (11) 

i|2 ,-, 



/ = 1 - \f'(zf 



■ip(«fti.or 

b(z;0;g 3 )l 2 
|t>fcft)l 2 



(15) 



is the inverse magnification of the individual image at the limit of 
the point source. In addition, the critical curves can be found from 
\f'(z)\ = 1 or equivalently solving f'(z) = -e~ 2 '* for z((f>) with 
being a real parameter (Witt 1990; Wambsganss et al. 1992; An & 
Evans 2006), and the caustics from their images under the mapping 
given by the lens equation. Since f'(z) is elliptic (hence, biperi- 
odic), there exists an infinite set of critical curves and caustics al- 
though they can all be recovered from the discrete grid translation 
of a 'unit' curve. 

The homogeneity relation of ^-function indicates that 



/ Cl)2\ 2 I Oil \ 

p(r,Q,K>)= — P{— 5 0,1) 



and from equation (A2) that 

e -2«i/3 / W2 \ 2 r /p™ i- 
! "'V3 



[-co 



(16) 



(17) 



Hence, the topography of the critical curves may be studied 
from the lines of constant values of \p(z; 1,0)| (square grid) and 
\p(z;0, 1)| (triangular and hexagonal grid). In particular, the criti- 
cal curves for the square (or triangular) lens lattice with the half- 
period oj are basically the curves of constant \p(z', 1,0) | = a> 2 /to^ 
[or \p(z; 0, 1)| = ll) 2 /^) 2 .] rescaled by a factor of w/ojo (or cj/to-i). On 
the other hand, those for the honeycomb lattice are found from the 
curves of constant \p{z\ 0, 1)| = a> 2 /( V3a>) scaled by V3w/ai 2 and 
rotated by -30°. Since p-functions are well-defined elliptic func- 
tions, for all three cases the critical curves are an infinite set of 
nonintersecting closed curves except when ll> assumes a particular 
value at which the curves are given by infinitely extended boundary 
lines that divide the lens plane. 

3.1 square lattice 

For the square grid case, if w < (o)o/V2) ~ 1.311 [i.e, > 
7r/(2wJ) ~ 0.457], each of the critical curves is centred at the centre 
of a square cell defined by the four adjacent lenses (i.e., kco + ipu> 
where k and p are odd integers). Any image in the region 'within' 
the critical curves has positive parity and vice versa. On the other 
hand, if a> > 2~ ,,2 ll> [k+ < 71/(2^)], the critical curves are cen- 
tred around each lens location. For this case, the images 'within' 
the critical curves now have negative parity. Finally, if a) = 2~ i,2 ll>q 
[k+ = 7t/(2Wp)], the critical curves are given by two infinite sets of 
parallel diagonal lines, and the whole lens plane is evenly divided 
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Figure 5. Critical curves (dotted lines) and caustic networks (solid lines) for 
lenses on a equilateral triangular lattice with cj = 1.25 [k* = 7t/(2 V3u 2 ) ~ 
0.580]. 



into checker-like tiled regions according to the parity of the images 
in them. 

For oj 4- 2~ ll2 a>o, the lens equation maps each connected and 
closed 'unit' critical curve to a caustic that is also connected and 
closed. Fig. 2 shows the critical curves and the caustic network 
for oj = 1.4. Like the critical curves, the whole caustic network 
is composed of the lattice arrangement of 'unit' caustic curves, 
but the unit curve exhibits self-crossing in contrast to the critical 
curves. The shape of the unit caustic curve, which is roughly sim- 
ilar to a variant of (irregular) octagram (i.e., {8/3)-star-polygon 3 ), 
is actually generic for the case that oj > 2~ [I2 ojq. The 'centre' of 
each unit caustic curve forms a similar square grid to the lens lat- 
tice but the side length of the grid on the source plane is given by 
2w(l - k+). As oj gets larger (or equivalently k+ gets smaller), the 
size of the caustics shrinks and they asymptotically reduce to de- 
generate points. For oj < 2~ i/2 oj , the caustic network can still be 
understood as the lattice arrangement of 'unit' caustic curves, the 
shape of which is simply described as being a diagonally stretched 
square. Again, the separation between the centres of the unit caus- 
tics are given by 2o> 1 1 - <c*|. However, this is smaller than the 'size' 
of each unit caustic if k+ ~ 1 so that the network can exhibit over- 
lapping among the neighboring caustics, which leads to a rather 
complex network (albeit regular thanks to the symmetry of the sys- 
tem). The critical curves and caustic network for co = 1.2 are shown 
in Fig. 3. The transition of the caustic network over oj = 2~ [/2 oj 
can be understood as the contact between (the vertices of) adjacent 
quasi-octagram caustics leading to their connection. The further 
evolution of the caustic topography for oj < 2~ ,,2 ojo is presented 
in Fig. 4. Since the separation between the unit caustics reduces to 
nil as K t — > 1, the network is the densest around K t « 1. On the 
other hand, as k+ increases past unity, the separations between unit 



A [p/^J-star-polygon, with p, q positive integers, is a figure formed by connecting with straight 
lines every (7-th point out of regularly spaced p points lying on a circumference. See http : 
//mathworld .wolfram. com/StarPolygon . html. 
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Figure 6. Same as Fig. 5 except for oj = 1.15 [at* = 7t/(2 V3w 2 ) » 0.685]. 

curves 2<jj(a:* - 1) increase while their sizes continually decrease. 
Consequently, the network eventually reduces an array of separate 
tiles of the unit caustics, which exhibit no overlap between neigh- 
boring caustics. The critical value at which the sizes of the unit 
caustics coincide with the separations between them is found ap- 
proximately to be oj x 0.748 Or* ~ 1.405). 

One physical interpretation of the caustics in gravitational 
lensing is that they are boundaries between regions in the source 
plane that produce different number of images. In addition, it is 
also known that a pair of images of opposite parity appears or dis- 
appears whenever the source crosses the caustics, and that the num- 
ber of positive parity images, if the source lies outside any caus- 
tics, is one or nil depending on the characteristics of the system. 
Hence, the maximum number of positive parity images can be de- 
termined from examination of the caustics network. For example, 
the quasi-octagram caustics for oj > 2~ xl2 oj( t allow at most four 
positive parity images in the region around the centres of caustics. 
It is easy to convince oneself that the minimum number of posi- 
tive parity images for the corresponding scenario is one when the 
source lies outside any caustic. As for the case oj < 2~ 1/2 oj , the 
examination of the caustic network leads us to conclude that, while 
the maximum number of the positive images approaches infinity at 
= 1 (i.e., oj = 7t l/2 /2 ss 0.886), it is greater than four only if 
0.803 < oj < 1.039. In other words, the number of positive parity 
images is bounded by four provided that oj > 1.039 (k* < 0.727) or 
that oj < 0.803 (a:* > 1.217). The source positions with no positive 
parity image start to exist if oj < 0.776 or equivalently > 1.039. 
Finally, if oj < 0.748 (k+ > 1.405), there is no overlap between 
neighboring caustics, and thus the number of positive parity images 
is one if the source lies in the caustic, and nil if otherwise. 

3.2 equilateral triangular lattice 

The lenses on an equilateral triangular lattice produce the critical 
curves and the caustic network in a basically consistent pattern as 
the square lattice lens although many details are quite different. For 
triangular lattices, the unit critical curves forw < (oj 2 /y/2) a 1.214 
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Figure 7. Critical curves (dotted lines) and caustic networks (solid lines) for 
lenses on a honeycomb lattice with a> = 1.45 [k+ = 7t/(3 V3<a< 2 ) ~ 0.286]. 



Figure 8. Same as Fig. 7 except for cj = 1.3 [** = tt/(3 V3w 2 ) ss 0.358]. 



[k* > 7t/(2 ' 3 3 1/2 <u|) ~ 0.615] are located in each of the equilateral 
triangular cell defined by three adjacent lenses. On the other hand, 
those for to > 2~ l,i t02 [f* > 7i/(2 ,l3 3 >l2 co 2 l )] are centred at each 
lens position. Note that, in contrast to the square lattice, the number 
of positions with zero shear is twice that of the poles (i.e., the lens 
positions) so that there is two-to-one correspondence between the 
unit critical curves of the former and that of the latter. At the crit- 
ical value that to = 2~ 1/3 toi2 [*<* = 7t/(2 1/3 3 1/2 wf)], the lens plane 
is divided into two classes by the critical curves: triangular regions 
around the null shear positions within which the images have posi- 
tive parity and hexagonal regions around the lenses corresponding 
to the images of negative parity. 

For to > 2~ 1/3 o)2, the shape of the unit caustics may be de- 
scribed as snowfiake-like or an irregular quasi-dodecagram (j 12/5}- 
star-polygon). Each unit caustic is arranged in the same equilat- 
eral triangular grid with side length of 2to(l — k+) to form the 
whole caustic network. Fig. 5 shows an example of this, the caus- 
tic network for to = 1.25. The maximum number of positive par- 
ity images associated with these unit caustics is found to be six. 
If to < 2~ l/3 a>2, the unit caustics are given by stretched triangu- 
lar cells, and the whole network is given by tiling the source plane 
with these cells, which may overlap with one another. Here, the 
centres of each cell actually form a hexagonal (honeycomb) grid 
with side length of 2a> |1 - / V3, which is in fact the mapped 
image of the dual grid to the lens lattice in the lens plane, and 
the adjacent cells alternate their orientation with their symmetry 
maintained. An example of these caustics networks is presented in 
Fig. 6 for ai = 1.15. Analogous to the case of the square lattice, 
dense overlapping of unit caustics leads to very complex network 
around ~ 1. It is found that that the critical value that corre- 
sponds to no overlap between adjacent caustics is about to « 0.791 
(k+ ss 1.449). The triangular lens lattice with side length less than 
this value (hence a denser lattice) produces triangular unit caustics 
arranged on a honeycomb grid with no overlap, and the number of 
positive parity images is either one or nil depending on the position 
of the source relative to the caustics. 



3.3 regular hexagonal lattice 

In qualitative terms, the critical curves of the honeycomb lattice are 
basically the 'dual' of that of the triangular case. The critical curves 
for co < (2 2/3 3~ 1/2 w 2 ) « 1.402 [** > tt/(3 1/2 2 4/3 w 2 ) « 0.308] are 
analogous to those for the triangular lattice with to > 2~ 1/3 W2 and 
located within each hexagonal cell centred at the positions of null 
shear. On the other hand, if to > 2 2,3 3~ >,2 (±>2, the critical curves 
are like those for the triangular lattice with to < 2~ 1/i co 2 and cen- 
tred around each lens positions. The dual characteristic further in- 
dicates that there are twice more poles than the locations of null 
shear (which are in fact fourth-order zeros), and thus there is one- 
to-two correspondence between unit critical curves of the former 
to the latter (formally, the 'winding number' of the former is twice 
that of the latter). 

On the other hand, the caustics of a honeycomb lattice are 
qualitatively distinct from those of a triangular lattice. Fig. 7 shows 
an example for a sparse honeycomb lattice with to = 1.45, which 
is in fact archetypal for to > 2 2l3 3~ ll2 to 2 . The unit caustics are 
given by a six-sided closed self-intersecting curve, which may be 
labeled a bi-triangle, and they are arranged in a similar hexagonal 
pattern as the lattice of the lens with side length 2ai(l - /c*). We 
also note that the maximum number of positive parity images as- 
sociated with these caustics is three. For a dense honeycomb lens 
lattice (to < 2 2/3 3~ 1/2 w 2 ), the unit caustics are found to be in the 
shape of a stretched hexagon (an example of which for to = 1.3 is 
given in Fig. 8) and arranged in a triangular grid with side length of 
2 y3a>\l - k*\. They also exhibit trends of varying overlap similar 
to the square or triangular lattice, as k+ gets larger past the unity. We 
find that the critical value for no overlapping caustics for the hon- 
eycomb lattice is approximately given by to ss 0.622 (k* ss 1.565). 

Despite their variety, we also notice some common character- 
istics of the caustic networks of three kinds of regular lensing lat- 
tices. For example, with an optical depth that is greater than a cer- 
tain critical value, the unit curves are more or less in a similar shape 
as the lensing lattice except for the fact that they are distorted in a 
way stretching the vetices radially outward. They densely overlap 
with one another around k+ ~ 1 . However, as K+ gets larger past the 
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Figure 9. Top: mean number of positive parity image for lenses on a square 
lattice (solid line), a triangular lattice (short-dashed line), a hexagonal lattice 
(long-dashed line), and for the random star field (dotted line). Also plotted 
in dot-dashed lines is the mean total magnification (/i) = (Y-k+ Y 2 . Bottom: 
mean number normalized by the mean total magnification. The line types 
are the same. 



unity, their size shrink whereas the separations among them grow, 
and so the networks eventually reduce to a simple array of tiles. 
Likewise, towards the lower optical depths below the critical value, 
the networks are comprised of star-shaped unit curves whose exact 
characteristics are related to meeting pattern of the lensing lattice 
at the common vertex points. 



4 MEAN NUMBER OF POSITIVE PARITY IMAGES 

The area in the caustic accounting for its multiplicity can be un- 
derstood as the covering fraction of the source plane by the region 
with 'extra image' pairs once it is properly normalized. In general, 
however, the normalization is ill-defined for localized lenses. This 
difficulty is ameliorated with the introduction of the lattice of point 
masses since the whole plane is now evenly divided into identical 
cells. With the regular lensing lattice, the area under the unit caustic 
can be properly normalized with respect to the area of the mapped 
image (on the source plane) of the unit cell comprising the lens- 
ing lattice. The result, if the multiplicity is properly accounted for, 
is actually the mean number of extra image pairs, which has been 
calculated for the case of the random star field using techniques 
from the probability theory (Wambsganss et al. 1992; Granot et al. 
2003). 

For 'dense' lattices that produce nonintersecting (but possibly 
overlapping) unit caustics, the relevant area under the unit caus- 
tics is straightforward to calculate since, despite overlapping, they 
are simply reduced to a set of closed loops. Then, the area under 
one such loop normalized to the area of a unit cell in the source 
plane, *C(1 - /e*) 2 where stf^ is the area under the unit cell de- 
fined by the adjacent lens positions in the lens plane, is the mean 
number both of extra image pairs and of positive parity images for 
a given source position because the corresponding critical curves 
in the image plane completely enclose the region where any pos- 
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Figure 10. Distribution of lensing shear, y. The line types are the same as 
Fig. 9. 



itive images reside and thus there are no positive images formed 
'outside' of any caustics. On the other hand, for 'sparse' lattices 
that produce the network composed of a grid set of self-intersecting 
'star-shaped' curves, the self-crossing of the unit structures requires 
some consideration on the meaning of 'area' that properly accounts 
for the image multiplicity. We argue that the solution is rather sim- 
ple. Once the area calculation is reduced to the line integral along 
its boundary (see e.g., An 2005; An & Evans 2006) using the fun- 
damental theorem of multivariate calculus 4 , the corresponding line 
integral along the caustic actually results in the area counted with 
multiplicity. That is to say, they can be considered as the 'area- 
excess' or the 'over-covering factor' of the source positions that 
produce positive parity images. Unlike the previous case, the nor- 
malization cell in the lens plane is the unit cell of the dual of the 
lensing lattice, i.e., the cell defined by the adjacent null shear loca- 
tions (which is in fact the centre of each cell of the lensing lattice). 
The corresponding critical curves lie completely within this unit 
normalization cell and furthermore enclose the regions of the im- 
age plane where negative parity images reside. Subsequently, there 
is at least one positive parity image for any given source position 
- that is, the source 'outside' caustics forms one positive parity 
image. As a result, the mean number of extra image pairs is less 
than the number of positive parity images by one. We note that the 
mean number of positive parity images actually varies continuously 
over the critical value dividing the dense lattice from the sparse one 
while the mean number of extra image pairs jumps by one. In fact, 



4 

The theorem is usually known as Green's theorem in elementary multivariate calculus when refer- 
ring to the relation between two-dimensional integral over a domain in two-dimensional space and 
one-dimensional integral along the boundary of the domain. This is also a special case of so-called 
(generalized) Stoke's theorem, the name under which the theorem is usually known as. 



© 2007 RAS, MNRAS 000, 1-19 



Lensing by regular point-mass lattices 9 



designating a certain image pair as an extra involves some degrees 
of arbitrariness whereas counting the number of positive parity im- 
ages is well-defined. For these reasons, henceforth, we shall discuss 
the mean number of positive parity images exclusively. 

Fig. 9 shows the resulting mean number of positive parity im- 
ages as a function of k+ . Comparison to the case of the random star 
field (with zero external shear) (Wambsganss et al. 1992) reveals 
some significant differences between the regular lensing lattices 
and the random field, particularly in the behaviour toward — > oo. 
While the mean number falls off as k~ 4 for the random star field and 
triangular lensing lattices (albeit the normalization for the random 
field being larger by a factor of ~4), its behaviours for the lenses on 
square lattices and honeycomb lattices are characterized by slower 



fall-offs given by 



and 



respectively. At first, this is 



somewhat counterintuitive. That is, in the random star field, the 
distance s to the closest neighboring lens is distributed according 
to 



&>(s)ds = 2ncrse- ncrs ds 



(18) 



where <x = K*/n is the mean surface number density of the stars. 
The mean and the variance are given by u> = 1/(2 y/cr) and [(4/n) - 
l]a) 2 = [1 - (7t/4)]k~'. In other words, as the density of stellar field 
increases, the random fluctuations around the mean decrease and 
one would expect that the random stellar field may approach to a 
regular field. 

The actual result can be understood using the distribution of 
the (internal) lensing shear y due to stars experienced by the given 
line of sight towards a lensed image. For a field of point masses, 
once this distribution is known, the mean number of positive parity 
images is recovered through (Wambsganss et al. 1992) 



(N + ) = </i> 



Jo 



y 2 )0»(y)dy. 



(19) 



With lenses on a regular lattice, we have shown that the shear at the 
given image plane location is given by the Weierstrass elliptic func- 
tion, and hence, the distribution on the whole image plane reduces 
to that on the unit cell in the image plane. Then, the cumulative 
distribution of shear is found by the fraction of area in the unit cell 
where the shear is smaller than the given value, and consequently, 
its derivative gives the differential distribution of the shear. Further- 
more, the homogeneity relation of ^-function implies the similarity 
relation between the distribution corresponding to different a>. For 
example, with a square lensing lattice, we have that 

£ ^(r'k.W = ^-*£,[|p(z;*2,0)| < y] 



■f 

Jo 



wo 



>(-z;i,o) 



7t y 



^(r'ko)dy' 



1 _d 

d r' 
i dg d 



(20) 



(21) 



*0 



where k = n/(2w ) 2 and ^Kj[\p(z; gi,0)\ < y] is the area of the 
unit cell where \p(z; g2,0)\ < y for a square lensing lattice with side 
length of w (see eq. [5] for g 2 as a function of u> and k*). In fact, any 



value of K-t instead of kq can be chosen as the normalization and the 
resulting similarity relation is valid for all three regular lattices. 

The resulting distributions for all three regular lattices are 
shown in Fig. 10 together with the same distribution for the random 
star field (Nityananda & Ostriker 1984; Schneider 1987; Wambs- 
ganss et al. 1992), which is given by 

& m (V\K*) = , 2 K *\, V2 = \- ^ T-375-- (22) 

(^+ r 2 )3/ 2 K\ [1 + (y/tfj 2 ] ' 

Unlike the random field case, we find that there is nonzero proba- 
bility that the line of sight experiences no net shear for the square 
lattice lens. By contrast, the same distribution linearly tends to zero 
as y — > for the triangular lattice (and also the random field) and it 
diverges as ~ y~ l/1 for the honeycomb lattice. These behaviours are 
related to the fact that the zeros of p(z; g 2 , 0), p(z; 0, #3) and b(z; g 3 ) 
are second, first, and fourth-order, respectively; that is to say, their 
Taylor-series expansions at their zeroes are given by 

p(z;g 2 ,0) - ( Z - z ) 2 + ffdz- z \ 6 ) 

p(z; 0, g 3 ) - igf (z - zo) + m\z - zol 4 ) (23) 

where zo is a zero of the corresponding function (note that zero off) 
is at zo = 0). For |y| <sc 1, the discussion in the preceding paragraph 
indicates that ^ 3 g (y')dy' =s n \z - zol 2 where zo is the null 
shear position and z traces the locations at which the shear is given 
by 7. By noticing the shear is given by | p(z; g 2 , 0)|, \p(z\ 0, g 3 )\ and 
Ibfc gi)\ for square, triangular, and honeycomb lattices, we can then 
find the leading-order behaviours for the shear distributions as 7 — > 
from the above Taylor expansions 5 ; 



&(rk*) - — 



1/4 

1 2 11 7t < ' 



^o(rk*) - 



1 2 4 7l"/ 2 

k* 3 3 / 4 r? 



1/3 



1/2 



0.1705 



(24) 



H-) - 



which agree with the numerical results shown in Fig. 10 (note that 
k+ = 7t for Fig. 10). In other words, the higher-than-linear-order 
(i.e., degenerate) zeros of the shear lead to nonzero finite or even 
divergent probability for the zero net shear. Physically speaking, 
accidental canceling of the net shear in the random star field is 
much less likely to occur than in the highly symmetric arrange- 
ment of the lenses, for which there may be significant chances that 
the image is located at an exactly balanced position where the net 
shear is null. This is basically the reason for slower fall-offs of the 
mean number of positive parity images in dense square or honey- 
comb lattices. From equation (21), if the shear distribution behaves 
like ^(jIa:*) - S^Hy/**)" as y/K t — > where S is a constant, 
a = (2 In) - 1 > -1 is the asymptotic power index for the shear dis- 
tribution and n is the order of the zeros of the shear, then equation 
(19) as k* — > 00 reduces to 



1 r 1 

+> = 71 3 d 7 (l- r 2 )^(7kJ 

(1 " K*) 2 JO 

f dy(l-yV. 
Jo 



(25) 



It is actually easier to work out from the series for the inverse functions, which are in fact hyper- 
geometric function (see also Appendix). 
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Figure 11. Locations of lensed images with a square lensing lattice with at = 1.35 (*■* » 0.431). The locations of the source with respect to the caustic 
networks are indicated in the top right panel. In the remaining three panels, the locations of the individual lensed images are marked by open circles (positive 
parity images) or closed dots (negative parity images). Also drawn in dotted lines are the critical curves and the locations of the lenses are marked by crosses. 



Since the last integral is a finite constant for a > -1, we find that 
(N+) ~ k* <1 * +3> as K-t — > oo, which is consistent with the results 
found earlier (Fig. 9). 

In fact, the shear distributions for regular lensing lattices are 
quite distinct from that of the random field, and it is somewhat re- 
markable that the behaviours of (N + ) are generally similar for all 
cases at least for small k+. This is related to the fact that (N+) is ba- 
sically an integration (or convolution) of the shear distribution, and 
therefore any 'roughness' in the shear distribution is 'smoothed' 
out in (N+). In particular, we note that the similarity relation of the 
shear distribution is essentially the result of the scale invariance and 
so should be generic. Then, &(y\K*) = K~ l< P(y/ic*) where ViyjKt,) 
is a function of y/*r*, and consequently we have 



(N + ) 



(1 



1 r l /K+ 



df(l - K 2 k f)<P(t)- 



(26) 



Furthermore, since the shear for the lines of sight that pass suffi- 
ciently close to any of the lenses is simply given by y a ct 2 where 
d is the separation between the given line of sight and that to the 
lens, it is straightforward to establish that the generic asymptotic 



behaviour of the shear distribution as y/K* — > oo to be £?{y\Ki,) — 
K+ly 1 [i.e., P(t) ^ r 2 as t — > oo]. Suppose that power-series expan- 
sion of P(t) at t = oo is given by P(t) - r 2 (l + 2,?=i c„r p ") (and it 
is uniformly convergent in an interval containing / = oo). Then, for 
k* — > 0, we find that (assuming p t l/n for any positive integer ri); 



Jo jri pn + 1 

Jo k* j~t P n ~ 1 



(27) 



2c„ 



pn+\ 



=! 1 - 2K+ + PnKl + > ■ 

(II) * * Z-i (p„ - l)( pn + 1) 



where P is a finite constant that weakly depends on the global 
(that is, far from / = oo) behaviour of V(f), and thus (N + )/(p) ^ 
1-2k* + ff(K 2 t ) and (N+) = 1 + 6\k\), regardless of the higher-order 
behaviour of V(t), provided that p > 1. In fact, we have < P rdn (t) = 
r(l + t 2 )- 3 ' 2 = r 2 (l + r 2 )- 312 so that p = 2 for P tw (f), and also 
the power-series expansions of ^-functions at their poles and zeros 
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(or equivalently the hypergeometric function form of their inverse 
functions) imply that the power-series expansions of P(f) at t = oo 
for regular lensing lattices are given in the above form with p = 2, 
p = 3, and p = 3/2 for the square, triangular, and honeycomb 
lattices, respectively. Therefore, all four cases exhibit the common 
asymptotic behaviour of (N+) for <k 1. 



4.1 image locations 

While the overall characteristics of the lensing situation can be un- 
derstood simply by studying the critical curves and caustics, one 
still needs to invert the lens equation to find the image locations 
for a given source position. With transcendental functions involved, 
this must be done numerically, but some analytical insights can ease 
the task. 

Since the real axis of the system for the lens equations (6), (8), 
and (11) is chosen such that the lensing lattice exhibits reflection 
symmetry with respect to this axis, the deflection functions for all 
cases should behave f(z) = f(z), which is indeed confirmed by the 
symmetry of the Weierstrass functions. Then, the lens equations 



2 4 -4 -2 2 

Figure 12. Same as Fig. 11 except a) = 0.8 (K t a 1.23). 

(6), (8), and (11) can also be written to be 



n=z-m-, m 



<T(z;0,g 3 ) a 
fcfcgs) o 



(28) 



Applying (two-dimensional) Newton-Raphson method to equation 
(28) indicates that the solutions are given by the (n — » oo)-limits of 
the sequences given by the recursion relation 



1 + /fa) -lZ n -fj- /fa) + Z„f'(Z„)]f'(Z„) 

l-/'fa)/'fa) 



(29) 



Alternatively, we can eliminate z from equation (28) by means of 
their complex conjugate (e.g., Witt & Mao 1995; An & Evans 2006) 



(30) 



and consequently, the image positions corresponding to the source 
position of 77 are the fixed points of the 'conformal' mapping given 
by 



Hz) = n + m + f(z)] ; h'(z) = f'[fj + f(z)] f'(z) 



(31) 
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Figure 13. Same as Fig. 11 except a) = 0.7 (k+ ss 1.60). Shown in the top left panel, there is no positive parity images for the source lying 'outside' the 
caustics. 



where 



f'(z) : 



-p(z;g 2 ,0) □ 
-p(z;Q,g}) a 
-b(z;g 3 ) o 



(32) 



In principle, all image positions can again be found applying the 
(complex) Newton-Raphson method to zeros of z - h(z) (which in- 
clude all image positions and possibly some spurious solutions) un- 
less the image is right on the critical curve. Specifically, the zeros 
can be located from the limits of the sequences given by 



Zn-h(Zn) Z„h' (Zn) ~ KZn) 



1 - h'(z„) 



h'(z„) - 1 



(33) 



as n — > oo. While equation (33) is basically the same as equa- 
tion (29) except for the fact that z„ in equation (33) is replaced by 
fj + f(z„) in equation (29), the latter relation in general is slightly 
better-behaved and also converges faster than the former. The it- 
eration (with infinite precision) given by equation (33) necessarily 
converges to a zero starting from anywhere in the complex plane 
except for a set of points with zero measure ('Julia set'), and fur- 
thermore any solution that is not on the critical curve has a basin 



of attraction that contains an open neighborhood of the solution. 
However, where the iteration converges to for a given starting point 
is nontrivial to figure out a priori. 

In fact, if we want to locate all of positive parity images, there 
actually exists an alternative route which is much more economi- 
cal. Application of the results from the complex dynamics (see also 
Appendix of An & Evans 2006) implies that, for any solution zo of 
equation (30) that is also |/'(zo)l < 1, there exists a point z c such 
that f'(zc) = and lim, M0O h"(z c ) = z - Here, h\z) = h(z) and 
h"(z) = h(h"~ l (z)), that is, h"(z) is the n-times iterations of the map- 
ping h starting from z. In other words, all the positive parity image 
positions are the limit points of the iterations of the mapping h(z) 
given by equation (31) starting from a zero of f'(z) (eq. [32]). Con- 
versely, the iterations of the mapping h(z) starting from any zero 
of f'(z) necessarily converge to one of the fixed points of h(z), the 
set of which contains all the positive parity image positions. Since 
the number of positive parity images can be determined by exam- 
ining the caustic network, one only needs, in principle, to perform 
the iterations finitely many times until all images are located (un- 
less k+ = 1). In addition, all the zeros of f'(z) for regular lensing 
lattices occur at the centre of each unit cell defined by the adjacent 
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Figure 14. Same as Fig. 1 1 except for an equilateral triangular lensing lattice with ai = 1.15 (a:* as 0.685). 



poles (here, the poles are simply the lens positions). In summary, 
all the positive parity image positions for a given source position 
can be located by the iterations of h(z) starting from the centres of 
each unit cell (corresponding to the null shear positions) until the 
set of distinct limit points contains the number of image positions 
of positive parity that is predicted by the source position and the 
caustic network. 

As for negative parity images, although there are infinitely 
many of them, we argue that most of them are insignificant. First, 
the lens equations indicate that \z — r}\ = |/(z)| = |/(z)|, which im- 
plies that, in order to have an image that is arbitrarily far from the 
source, the image position should be rather close to the lens, that is, 
the poles of /(z), at which /(z) is divergent. If we consider the lens 
mapping from the point z = 8z + fl m „ where Cl m „ is one of the lens 
locations, then the lens equation (28) reduces to 



7] = Sz + fi,„„ - f(Sz + n,„„) = Sz + (1 - **)fi„,„ - f(Sz) 



(34) 



where Sz = Sz- Here, we use the quasiperiodic relation of f(z) 
[which is inherited from that of ( (z)] and some additional proper- 
ties of the Weierstrass zeta function. Note also that k+ = nl (2cm) 1 , 
K t = n/(2 VSw 2 ), and k* = 7t/(3 V5w 2 ) for square, triangular, and 
honeycomb lattices with side length of w. Assuming |&| <k u, 



equation (34) can be solved approximately from 

(1 - - T] = f(6z) -5z^--6z + M 

Sz 



(35) 



where the remainder term (M) is the order of (6zY with j = 3 
(square), 7 = 5 (triangle), or j = 2 (hexagon). Provided that |Jz?| » 
oj~ [ where Jzf = |_Sf| e"* L = (1 - k^Q.,,,,, - r], we have the image at 
z = fi„,„ + Sz where 



and its signed magnification is found to be 
1 1 



(36) 



1 - \f'(z)\ 2 1 - \f'(Szf 
J_ . J_ / 1 \ 

|^| 4 + |^| 6 + V |ntoC/+5,8) j' 



(37) 



Since Cl m „ — fj = J^f/(1 - /c*) where /) = 77/(1 - ^*) with ^ 1, 
the condition that > L is sufficiently large translates to the 
lens positions (fS m „) that lie outside the circular region of radius 
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Figure 15. Critical curves and caustics for square lattice lenses with exter- 
nal shear (solid lines). Also plotted in dashed lines are the critical curves 
and caustics of the Chang-Refsdal lens (a single point-mass lens) with the 
same external shear. All four panels are for \y c \ = 0.25, which is oriented 
parallel to the side of the lattice cell, but varies; 0.1 (top left), 0.2 (top 
right), 0.3 (bottom left), and 0.35 (bottom right). 



Lo/(l - k*) centred at f}. 6 Then, the images near all of those lens 
positions are simply found by the approximation (eq. [36]). More 
importantly, the contribution to the total magnification from all of 
those images is approximately 



m,n m.n 1-^1 ° " 



2-no-RdR 



wci-**) |-£f| 
Ik+RAR 

ifCl-K*) 4 



(38) 



where the summation is over the lens positions (m,n) such that 
~ f]\ > L Q I |1 - K *\ and cr = K+jn is the number density of 
the lenses. Here the polar coordinate centred at r\ is used for the in- 
tegral and soR = \Q.„,„ - f)\. The result indicates that, unless k+ = 1, 
the contribution from infinitely many negative parity images to the 
total magnification is finite 7 (note that the integrand is convex so 
that the sum is bounded by the integral) and with sufficiently large 
Lo, they are negligible. 

In Figs. 1 1-14, we show the image locations for some specific 
cases of lattice lenses. Here, all positive parity images are located 
through the iteration scheme outlined earlier. Most of negative par- 
ity images are found via the Newton-Raphson method (eq. [33]) 
with the starting locations given by the approximate solutions in 
equation (36). Additional locations of negative parity images are 
also found from the iteration given by equation (33) with starting 
points given by a grid of positions near z = 77/(1 - k+). 



Note that, as K* — * 1, the radius of the region within which the approximation is not valid also 
increases. 

The sum Yj'mn where the summation is over the lattice points except the points at which 

Jzf = 0, is 'absolutely convergent' for r > 2. For the rigorous proof, see the mathematical references 
on elliptic functions listed in Appendix. 



5 EXTERNAL SHEAR 

Next, we consider the effect of the large-scale potential. It is a usual 
practice that the external potential (fr Mt in equation (2) is approxi- 
mated by a quadratic function and consequently the correspond- 
ing deflection by a linear function (i.e., tensor) - that is, keeping 
only the linear terms in its Taylor-series expansion (Young 1981; 
Nityananda & Ostriker 1984; Subramanian et al. 1985). The result- 
ing linear function is symmetric, thanks to the mixed coordinate 
differential being commutative. It is also customary that the linear 
function is decomposed into a scalar dilation and a traceless ten- 
sor. The former, usually referred to as a convergence, is directly re- 
lated to the surface mass density (which is not in the form of point 
masses) whereas the latter, known as an external shear, physically 
models the tidal effects from the external mass distribution. 

In complex number notation, then the lens equation general- 
izes to 



r] = (1-k c )z-j c z- f(z); 



V-Kcf- [-/'(z)]-rJ (39) 



where k c , which must be positive real, is the convergence due to 
the local surface mass density that is the source of the large scale 
potential (i.e., 'continuous mass density') and y c , which can be any 
complex number, is the shear that models the tidal effect. The ef- 
fect of the convergence term is uniform focusing, which basically 
results in a scaling difference between the image and source planes. 
It is usually incorporated into the analysis by the redefinition of the 
variables (Paczyiiski 1986). For example, with the rescaled posi- 
tion variables for the image plane w = |e| I/2 z and the source plane 
£ = sgn(e) |eP I/2 rj where e = 1 - k c , equation (39) for the square 



lattice with side length of to [i.e., f(z) 
reduces to 



ftegz.O) = f(*|<u,ku)] 



f = w - gw- sgn(e) \e\ L ' 1 oj,\\e\ 1 ' 1 oj) ; 

7 =1-Uw| |e| 1/2 ^,ik| I/2 ^)-sgn(e)d 2 



(40) 



a-K C ) 



where <; = y c /e is the reduced shear. We note that the deflection 
function is basically given by that for the square lattice with side 
length of |e| 1/2 td (i.e., corresponding to the reduced optical depth 
x = k+Ie). In general, if e > (i.e., < K c < 1), the qualitative 
lensing behaviour is basically same as the case that k c = except 
the additional scaling difference between the image and the source 
planes, which leads to the use of the 'reduced' values of the external 
shear and the optical depth (Paczyiiski 1986; Granot et al. 2003; An 
& Evans 2006). For the following discussion, we only consider the 
lens equation (39) for the case k c = 0. 

The critical curves and the caustic networks of the lattice lens 
with external shear can be found in the same way as the case with- 
out external shear. However, the shear introduces two further de- 
grees of freedom in the parameter space to be explored, which 
makes the task quite extensive and difficult to generalize. Never- 
theless, we can draw a rough generalization such that the patterns 
reduce to the case of the Chang-Refsdal lens (Chang & Refsdal 
1979, 1984; An & Evans 2006) if \y c \ » k* and to the simple lat- 
tices with no (or negligible) shear if \y c \ <K /e*. Figs. 15-17 show 
some examples of the critical curves and the caustic networks for 
the square lattices with moderate values of \y c \ and to illustrate 
the effects of varying k+ at fixed y c or varying y c (both the magni- 
tude and the orientation) at fixed K t . The situation becomes quite 
complicated if \y c \ « k* ~ 1. As in the case of lattices with zero 
shear, with /c* ~ 1, there are overlaps among neighboring caustics. 
Furthermore, with nonzero shear, neighboring caustics can be con- 
nected with one another, and the resulting dense networks of caus- 
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Figure 16. Same as Fig. 15, but with \y c \ = K t = 0.25. Four panels show 
different orientations of the external shear: 0° (top left), 22? 5 (top right), 
30° (bottom left), and 45° (bottom right). Note that the figures are rotated 
so that the direction of the external shear is along the horizontal direction. 

tics are even more complicated to analyze properly. More detailed 
analysis of these will be attempted in a subsequent work. 

Instead, here we speculate on some possible connections of 
lattice lensing to more realistic lensing scenarios (particularly, that 
of the random star field). One of the most significant differences 
of lattice lensing from that by a random field is the fact that the 
shear experienced at any given lens position due to the remaining 
lenses cancels out to be null. Adding external shear only alters the 
situation minimally, that is, the total shear at any lens position is 
shifted to the value of external shear but every lens still experiences 
the identical shear. By contrast, in the random star field, the shear 
experienced by a given star is distributed according to the same 
distribution as that of the shear experienced by the random lines 
of sight. However, we note that the latter distribution has been cal- 
culated (Nityananda & Ostriker 1984; Schneider 1987; Granot et 
al. 2003) using the straightforward application of the characteristic 
function in the probability theory and the random-walk process. If 
lattice lensing were to model the effect of the shear 'variations' in 
the random field reasonably well, each individual lens in the ran- 
dom field might be approximated as a lattice lens with the 'external 
shear' being equal to the actual total shear experienced by the given 
lens. Then, the lensing behaviour of the whole of the random star 
field could be understood as the sort of 'convolution' of the lattice 
lens (of the corresponding optical depth) with the external shear 
distributed as the distribution of the random shear in the random 
star field. If this indeed be the case, one could study more realistic 
scenarios of the microlensing effects due to the stars in the lensing 
galaxy in an alternative way to the traditional inverse ray tracing 
approach. 

At this point, however, this remains a speculative possibility. 
As it has been noted widely, the effects of many lenses combine 
in a highly nonlinear fashion and the role of the individual lenses 
even in the simplest phenomenology is nontrivial (see e.g., Gil- 
Merino & Lewis 2005). While one could expect that the approach 
of decomposing the star field into lattices outlined in the previous 




Figure 17. Same as Fig. 15, but with a fixed value of /c* = 0.25 and varying 
\y c \: 0.3 (top left), 0.1 (top right), 0.05 (bottom left), and 0.01 (bottom right). 
For all four, the orientations of the shear are parallel to the side of the lattice 
cell. 

paragraph might fare better than the simplistic approach of seeing 
the star field as the combination of individual stars (since it takes 
care of some of the nonlinear effects), no compelling argument can 
yet be made that this is the case. 



6 CONCLUSION 

Microlensing at high optical depth is not only of a theoretical inter- 
est but also of an importance in understanding a number of lensed 
systems (e.g, Keeton et al. 2006; Morgan et al. 2006). However, 
at present time, the only way to get any detailed quantitative grip 
on such problem is via Monte-Carlo simulations of each individual 
point corresponding specific values of the optical depth and the ex- 
ternal shear (e.g., Granot et al. 2003; Schechter et al. 2004), which 
is rather expensive in its requirement for resources. We argue that a 
grid of point masses discussed in this paper, while artificial and ab- 
stract so that it may appear to be entirely divorced from any phys- 
ically realistic problem, can provide us with some insight, albeit 
indirect, into what is at work in microlensing at high optical depth, 
in particular for the numbers and magnifications of microimages. 
Together with with various analytical methods and probabilistic ap- 
proaches to the problem, this help us build pictures of phenomena 
complement to more traditional routes with comparatively inexpen- 
sive means. 
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APPENDIX A: WEIERSTRASS ELLIPTIC FUNCTIONS 

Here, we summarize some basic concepts and properties of ellip- 
tic functions and the Weierstrass functions that are related to the 
expositions given in the main text. For a general but concise intro- 
duction to the subject, we refer to online resources such as Wol- 
fram Mathworld 8 or standard references on special functions such 
as Abramowitz & Stegun (1972). For more detailed and rigorous 
(and modern) treatments, we advise our reader to look upon ad- 
vanced texts on the complex analysis (e.g., Ahlfors 1979; Lang 
1999) or the analytic number theory (e.g., Apostol 1997), particu- 
larly those on the theory of the elliptic curves and/or modular func- 
tions/forms (e.g., Silverman 1986; Lang 1987; Koblitz 1993). In 
addition, Large collections of formulae involving the Weierstrass 
functions are also available online (e.g., the Wolfram functions 
site"). 



Elliptic functions are the class of complex-valued functions 
of a complex variable that are meromorphic (i.e., analytic every- 
where in the complex plane except isolated poles) and biperiodic 
(or 'doubly-periodic'). The biperiodicity means that there exist 
two complex numbers a>[ and o> 3 such that a>i/o) 3 is not real and 
f(z + 2mo) l + 2no) 3 ) = f(z) for all pairs of integers m and n and any 
complex number z. Here, the two complex numbers oil and o) 3 are 
referred to as the half-periods of the elliptic function if they are a 
pair of two 'smallest' such numbers. 10 We note that the biperiodic- 
ity further implies that, once the behaviour of the function is known 
in the domain that is the parallelogram 'cell' whose vertices are 0, 
2a>i, 2w 2 = 2a>[ + 2o) 3 and 2o) 3 (or equivalently oil - o) 3 ,o) t +0)3, 
-a)i + w 3 , -oil - o) 3 ), it is completely specified everywhere in the 
complex plane. The cell is sometime referred to as the fundamen- 
tal period parallelogram (FPP) of the elliptic function. According 
to the general theory of elliptic functions, any elliptic function can 
be expressed using a set of standard functions, which are usually 
chosen to be one of two classes: Jacobi functions (sni, cnx, dnx 
and so on) or the Weierstrass p-function and its derivative. 

The Weierstrass p-function (or 'the' Weierstrass elliptic func- 
tion) is the elliptic function that has one second-order pole in its 
FPP. The standard definition locates a pole at z = whose principal 
part is exactly z~ 2 and constant part is nil (i.e., lim z ^ [p(z) - z~ 2 ] = 
0). The formal definition is given by 



P(z\oj 1 ,o) 3 ) = + ^ 

^ m,n=-co 
(m,n)*(0,0) 



1 



(z- 



1 

mil 



where Q.,„„ = 2mo>\ + 2no) 3 , and the summation is over all inte- 
ger grid points of (m, n) except (m,n) = (0,0). Here, the set of 
all complex numbers in the form of Cl mn is usually referred to as 
a lattice. Note that the function p(z) has second-order poles at ev- 
ery point of the lattice. In addition, p(z) is an even function, i.e., 
p(-z) = p(z)- Note that, for a given lattice, the choice of the pair 
of the half-periods is not unique - that is to say, the lattice forms 
a module over the ring of integers and the set of two periods (i.e., 
{2oj!, 2<^i 3 }) is the basis of this module. For this reason, the two 
secondary arguments of p(z) are usually replaced by the elliptic in- 
variants defined by g 2 = 60 YI mfl f2 m * and g 3 = 140 Yl mfl ®-~ m \ such 
that p(z\g2,g 3 ) = p(z\o)i,o) 3 ). Unlike a pair of the half-periods, the 
elliptic invariants are in fact uniquely specified by the given lattice, 
and vice versa. 

The elliptic invariants are directly related to the coefficients 
of the Laurent-series expansion of p-function at origin via simple 
algebraic relations, that is, the coefficients are given by the poly- 
nomials of g2 and g 3 . Particularly, we have p(z;g2,g 3 ) - z~ 2 + 
(g2/2Q)z 2 +(g 3 /28)z 4 +0(z 6 )- This series expansion further indicates 
that [p'(z;g 2 ,g 3 )] 2 - 4[p(z; g 2 , g 3 )] 3 - g 2 p(z; g 2 , ft) = ft + ^(z 2 )- 
This combined with the results from the general theory of ellip- 
tic functions implies that y(z) = p(z;g 2 ,g 3 ) is a particular solu- 
tion of the differential equation given by (y') 2 = 4y 3 - g 2 y - g 3 . 
Since this is a first-order differential equation that does not in- 
volve the independent variable, its general solution is given by 
y(z) = p{z + C;g 2 ,g 3 ) where C is an integration constant (which 
is an arbitrary complex constant). Note that, if the polynomial 
4y 3 _ g2)> ~ g3 can be factored with a perfect square factor (so that 
g\ = 27gj), the solution of the differential equation then can be ex- 
pressed in terms of 'elementary' functions such as trigonometric or 



http : //mathworld . wolfram. com/ 
http : //functions . wolfram. com/ 



'° The reader should be warned that, while we follow the convention of Abramowitz & Stegun 
(1972) in the paper, that is, w; being half-periods, many modern mathematical references adopt the 
different convention such that ojj and the related arguments represent the periods. 
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hyperbolic functions and so does the Weierstrass function. In par- 
ticular, if g 2 = 12e 2 > and g 3 = -8e 3 are both nonzero real [i.e., 
■ giy -g3 = My - ef(y + 2e)l then 

(3ecoth 2 [(3e) 1/2 z]-2e 
(-3e)cot 2 [(-3e) 1/2 z] + (-2e) 



4y 3 ■ 
p(z;12e 2 ,-8e 3 ): 



if e > 
if e < 



However, this case actually corresponds (at least, formally) loa>[ = 
roorw 3 = oo, and the resulting functions are not elliptic but simply 
(singly-)periodic [unless g 2 = g 3 = 0, for which p(z; 0, 0) = z~ 2 
and W] = oj 3 = oo so that it is in fact aperiodic]. 
The derivative of p-function 



p'(z;g 2 ,g3) = p\z\co u co 3 ) = —p{z\oJum) 



Yj Yj ( z -a 



)3 
n) 



is also an elliptic function with the same half-periods as p(z\aii, 103) 
but it has third-order poles at the lattice points and is an odd func- 
tion. By contrast, the antiderivative (or the indefinite integral) of 
p-function cannot be an elliptic function. According to the general 
theory of elliptic functions, the sum of complex residues in FPP of 
any elliptic function must vanish. However, if the antiderivative of 
p-function were an elliptic function, it would have a single simple 
pole in its FPP and thus the sum of its residues could not vanish. In 
fact, the Weierstrass ^-function", which is defined to be a particu- 
lar anntiderivative of (O-function; 



((z;g2,gi) = f(z|«i,« 3 ) = ~ + J" [~2 " pM^w >.<)|du 

= 7 + Z'(^ 



n 2 



) 



is a quasiperiodic function such that f(z + 2ma) l + 2noj 3 ) = f(z) + 
2m£(<jj 1 )+2n£(a> 3 ) for all pairs of integers in and n and any complex 
number z with oji and 0)3 being the half-periods of the correspond- 
ing p-function (secondary arguments suppressed for brevity). We 
also note that the Weierstrass ^-function is an odd function. 

In addition, the theory of the Weierstrass functions also in- 
troduces another auxiliary function referred to as the Weierstrass 
cr-function; 

<r(z;g 2 ,g3) = oto*,**) . z JT [(l - 0*v(± + ^r) 

(m,n)*(0,0) 

where the product is again over all integer pairs (m, n) except 
(m, n) = (0, 0). From this definition, it is straightforward to show 
that f (z) is the log-derivative of <x(z), that is, 

""'fe) d 1 / n r, n 

— — = —l nC r(z\g2,g3) = l;(z;g 2 ,g3)- 

tr(z) dz 

We note that <x(z) is also a quasiperiodic function although the re- 
lation involves the successive ratios between the functional values 
at the points separated by the 'period' , rather than the differences 
as in £ (z); 

<x(z + 2mtL) l + 2na>3) 



<r(z) 



■ (-If exp {2[m(((jj) + n{(a>)](z + ma> 1 + noj 3 )} 



where p = m + n + mn and m and n are any pair of integers. The 
function <x(z) is in fact an entire function, that is, it is analytic every- 
where in the complex plane without any pole or singularity (except 
the essential singularity at infinity). In particular, it has zeros at the 



^ ^ The Weierstrass ^-function should not be confused with several other mathematical functions 
conventionally denoted by the same symbols, the most important among which is the Riemann f- 
function and its generalization, the Hurwitz ( -function. 



origin and the lattice points, that is, cr(0) = cr(fl,„„) = 0, which is 
easily verifiable from its definition. 

We also note that the set of Weierstrass functions satisfy the 
homogeneity relations, which can be easily derived from their re- 
spective definitions in terms of the infinite sum or product over the 
lattice points. In particular, 



cr(az\aLij],aaj 3 ) = acr{z\cO{, 0J3) ; 
£(az\aa)i,acL>3) = aT { ((z\oj 1,(1)3) ; 
p(az\au>i, au> 3 ) = a~ 2 <p(z|^i, 0)3) ; 



a-(az\ar A g2,ar 6 g3) = acr(z; g 2 , g 3 ) , 
t;(az\aT 4 g2,aT 6 g3) = fl"'f(z; g 2 ,g3) , 
p(az\aT 4 g2,ar 6 g3) = aT 2 p{z\g2,g3) , 



p'{az\a(D U acL>3) =a V(z|a>i, m ) ; p'(az;a 4 g 2 ,a 6 g 3 ) = a V(z; g2,g3) , 

g 2 (aoji,aaj 3 ) = a' A g 2 (a>i, co 3 ) , g 3 (aco u aa> 3 ) = a' 6 g 3 (a> u oj 3 ) , 
where a ^ is a complex constant. 



Al lemniscatic case 

If FPP is a square (i.e., co 3 = ia>i), then Q. mn = 2ma>\ + 
2nu>3 = 2(m + ni)a>i and, from the symmetry, we have Ylmn ^Knn = 
(2(x> l y 6 Y!mrA m + wl ) 6 = 0, and so that g 3 = 0. Then, after g 2 is 
reduced to the unity using the homogeneity relation, the behaviour 
of the corresponding Weierstrass elliptic function can be studied 
through that of p(z; 1,0), which is called the lemniscatic case. The 
corresponding half-periods of p(z\ 1,0) = p(z\d>Q, id>o) are given by 
a>o = r 2 /4 /(4 yfn) ss 1.854 and ioj . (The constant V2<x> ~ 2.622 is 
sometimes known as the lemniscate constant.) Some basic results 
known for the lemniscatic case are 



p(a> ; 1,0) = -: 
$5(zo;l,0) = 0; 



(iw ;l,0): 



f (wo; 1,0)= — , 
4o> 

f(z ;l,0) = ^-d-i), 

4w 

ai(o ;i,0) = -^-. 

4co 



Here, zo = (l+i)^o = V^e 711 ' 4 ^ is in fact the only zero of <p(z; 1,0) 
within its FPP (which is actually a degenerate second-order zero). 
Note that these zeros are located at the centres of the square cells 
defined by the poles. We further note that the behaviours of these 
functions within FPP exhibit high degrees of symmetry and their 
study can be reduced to the one-eighth of FPP, that is, the isosceles- 
right-triangular region defined by z = (the pole), co and zo (the 
zero). 

Although the standard reference on the numerical calcula- 
tion, Press et al. (1992), does not list the routine for the Weier- 
strass functions, their numerical evaluation is available in commer- 
cial softwares such as Mathematica® or Maple or subroutine li- 
braries such as IMSL 12 . However, for most purposes, they can 
be evaluated in a pretty straightforward manner without using any 
black-box routine (see e.g., Abramowitz & Stegun 1972). While 
there exist more involved algorithms that are valid for an arbi- 
trary pair of complex numbers (g2,g 3 ), the lemniscatic case of the 
Weierstrass function can be easily evaluated through its Laurent- 
series expansion at z = or Taylor-series expansion at z = w or 
z = Zo (after the reduction of the domain using biperiodicity and 
additional symmetry properties). The relevant series expansions up 
to the first several terms are available in Abramowitz & Stegun 
(1972). In particular, the Laurent series at z = is given in the 
form of p(z; 1,0) = z~ 2 [l + 2Xi(«*/ 2 0V*] (the first few coef- 
ficients a k 's of which are given in Table Al) so that it converges 
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Table Al. Laurent-Series Coefficients of p(z) at z = 0. 



k 


ak+i/ak 


h+i/h 


1 


1/3 


1/13 


2 


2/13 


1/19 


3 


5/(2- 17) = 5/34 


3/(5 ■ 13) = 3/65 


4 


2/(3 -5) = 2/15 


2 2 /(3-31) = 4/93 


5 


5/(3- 13) = 5/39 


(5 ■ 7 ■ 43)/(2 2 ■ 13 ■ 19 ■ 37) = 1505/36556 


6 


(2 ■ 3 2 )/(5 ■ 29) = 18/145 


(2 ■ 3 ■ 431)/(5 ■ 7 ■ 43 2 ) = 2586/64715 



For both cases, m = b\ = 1. 



relatively quickly for |z| < 1 - the radius of the convergence of the 
series is in fact 2ll> . Note that the coefficients for the Laurent series 
at z = can in fact be derived using a total recurrence relation. 

Since the lemniscatic case of p-function satisfies the differ- 
ential equation (y') 2 = 4y 3 - y, its inverse function is express- 
ible through elliptic integrals. However, the relevant integral in fact 
more easily reduces to an incomplete beta function 



:= ^ ;i '° )= -rj^4Mi4) 

H2'4'4'4<p 2 / 



r, 1,0) = 



J4f 

or a hypergeometric function, 

1 _ l\ 15 

l2'4'4'4<p 2 

where B x (a,b) is the incomplete Beta function and 2 Fi(a, b;c\x) 
is the Gaussian (2, l)-hypergeometric function, [c.f., B x {a,b) = 
(x a /a) 2 Fi(l —b, a; a + 1; x).] Here, the last expression may be used 
for a convergent hypergeometric power-series expansion for the in- 
verse p-function for \p\ » 2" 1 . We note that, because we are inter- 
ested in the case when p takes an arbitrary complex value, there is 
no particular advantage using the Legendre/Jacobi elliptic integrals 
or the incomplete Beta function, numerical or otherwise, compared 
to the hypergeometric series expression given above. Alternative 
convergent hypergeometric series for \p\ <K 2" 1 also exists; 

'1 1 5 



: = p-\p; 1,0) = z + 2ip 1 ' 2 2 F,(-, -; -;4p 2 ) 



where zo = v2e ra/4 Wo is again the zero of p(z; 1,0). The inverse 
function can be effectively evaluated almost everywhere using these 
hypergeometric series together with some further refinement meth- 
ods such as the (complex) Newton-Raphson algorithm if necessary. 

A2 equiharmonic case 

We note that FPP can be considered as a tile, the set of which cov- 
ers the whole (complex) plane and the locations of the poles are 
vertices of individual tiles. Then, the preceding lemniscatic case 
corresponds to the tile given by a square (i.e, a regular tetra-gon), 
which is one of the three possible regular tessellations of the plane. 
The elliptic function with FPP corresponding to the tile pattern of 
the remaining two regular tessellations can also be studied simi- 
larly. The half-periods for the more basic of the two, that is, the 
tiles given by equilateral triangles, are related to each other via 
u>3 = e ml3 a>i and then from the symmetry of the system, it fol- 
lows that g 2 = 602 m „'^™ = (note that, while the choice of the 
two half-periods is not unique, the elliptic invariants are fixed for 
the given arrangement of poles). As before, the remaining invariant 
g 3 can be reduced to the unity using the homogeneity relation, and 
the resulting Weierstrass elliptic function, p(z;0, 1) is referred to 
as the equiharmonic case. It is known that the real half-period of 



p(z;0, 1) = p(z\(0 2 ,e ni,3 (0 2 ) is w 2 = r 3 1/3 /(4n) ~ 1.530. Addition- 
ally, it can be derived that 

TTf.m/3 

p(e- m/ V ; 0, 1) = 2" 2/3 e 2m/3 ; ^ m/3 w 2 -0,1) 



2V3w 2 ' 



p(6; 2 ;0,l) = 2- 2/3 ; 
S5(e m/3 w 2 ;0,l) = 2- 2/ V 2m/3 : 
p(V3iw 2 ;0,l) = 2- 2/3 ; 



f(w 2 ;0,i) 



2V3w 2 



«e 7ti/3 ^ 2 ;0,l) = 



-m/3 



2V3( 



w 2 



«^;ai ) = -_. 



The zeros of p(z\ 0, 1) are found at the centres of each equilateral 
triangular cells so that there are actually two zeros within its FPP, 
and unlike the lemniscatic case, each zero is simple (i.e., the first- 
order zero). The particular zero within the cell defined by z = 0, 
z = 2w 2 , and z = 2e n,/3 oj 2 is located at zo = (2/ y/3)e n,l6 co 2 with 
<T(zo;0,l) = 7te-™ /6 /(3^ 2 ). 

Like the lemniscatic case, the equiharmonic case of the Weier- 
strass function can also be readily evaluated via the power-series 
expansions at z = 0, z = co 2 , or z = Zo- Their coefficients are 
again found in the references such as Abramowitz & Stegun (1972). 
In fact, the Laurent-series expansion at z = of the equihar- 
monic case (whose radius of convergence is 2oj 2 ) is in the form of 
p(z\0, 1) = z~ 2 [l + 2S=i(^/28*)z M ], which converges even faster 
than that of the lemniscatic case for |z| < 1. The first few coefficient 
b k 's are again given in Table Al. 

Analogous to the lemniscatic case, the inverse of p-function 
can be obtained through the integration of the differential equation 
(y') 2 = 4y 3 - 1. The resulting expression can again be written down 
in terms of elliptic integrals (which is rather complicated), incom- 
plete beta functions, or hypergeometric functions, 



j;0, 1) = - 



dy 



1 



B 



l/(4p 3 ) 



V47^T 2 2 / 3 • 3 

= _L F (l 11 -L) 

pi/ 2 2 l \2' 6' 6' 4sW' 
Similarly, the last expression can be used for the convergent hy- 
pergeometric power-series expansion of the inverse p-function for 
\<p\ » 2~ 2/3 . For the case when \p\ <k 2~ 2/3 , an alternative expres- 
sion involving a hypergeometric function can be derived as before; 

z = p-\p;0, 1) = zo - ip 2 Fi{^, p V). 

Here, zo = (2/ yj3)e ni,6 to 2 is again the zero of s o(z;0, 1). The 
sign in front of the hypergeometric function is actually depen- 
dent upon the choice of zo- For example, with an alternative zero 
Zo = -(21 V3)e 7tl/6 ai 2 , it should be positive. Again, these hyper- 
geometric series supplemented by the Newton-Raphson algorithm 
provide an efficient way to evaluate the inverse p-function for the 
equiharmonic case. 
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A2. 1 hexagonal grid 

Let us consider an elliptic function, 

Mr, ft) = P(z;0,g 3 ) - p( z ;0,-^). (Al) 

Note that every pole of p(z; 0, -g 3 /27) coincides with that of 
p(z\ 0, g 3 ) and they cancel each other leaving a fourth-order zero at 
the location. On the other hand, the remaining poles of p(z;0,g 3 ) 
are located at the vertices of regular hexagonal tiles. Hence, equa- 
tion (Al) corresponds to the remaining case of the elliptic function 
with their (second-order) poles forming vertex points of the regular 
tessellation of the space. 

While it is a straightforward calculation, equation (Al) is un- 
suitable for the accurate numerical evaluation near z = (and also 
near the other zeros of the functions) because of the canceling of 
two divergent terms leading to the zero. Furthermore, the inver- 
sion of equation (Al) is also rather nontrivial. This difficulty can 
be overcome as follows. We first note that 

W(z)f 



= [p'(z)] 2 (l- 2fi 
= (4[p(z)] 



27[p(z)] 

§)('- 



2g3 



27[p(z)] 



;) 2 = 4[q(z)] 3 



- g3 



where q(z) = p(z) + (# 3 /27)[p(z)r 2 and p(z) = p(z;0,-g 3 /27). In 
other words, q(z) is the solution of the differential equation (y') 2 = 
4y 3 —g 3 , whose general solution is in the form of y = p(z + C; 0, g 3 ) 
where C is the (complex) integration constant. Since the principal 
part of the Laurent- series expansion of q(z) at z = is z~ 2 , it is in 
fact the case that q(z) = f(z; 0, g 3 ), and consequently we have 
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*->=§[^o.-g)r=i(§)"i4€)"^') 

(A2) 

where u 6 = -1. Subsequently, we find that the inverse function of 
equation (Al) or equivalently that of equation (A2) is given by 

= 1,114 F (- - I 1)312 ) 

Z fe/27)^ 2 '{re' 6' 4(g 3 /27)i/ 2 J 



r 3 

1 1/3 



lFl l2'3'3' 



1,1/2 



4(g 3 /27) 



1/2 N 



where some of signs are chosen so that the principal value for pos- 
itive real values of b returns with being positive real (in particular, 
h > ^ z e [0, 2w)). 

On the other hand, if one wants to find the Taylor-series ex- 
pansion of b(z; g 3 ) at z = 0, it may be easier to use equation (Al) at 
z = directly; 



K= 1 L X ' J 



where b k is the Laurent-series coefficient of the equiharmonic case 
of ip-function (see Table Al). In particular, b x = 1 which is actually 
used for the second equality of the preceding equation. 



This paper has been typeset from a TjX/ MIfX file prepared by the 
author. 
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